Modeling catalytic combustion om methane using RMG-Cat¶

Based somewhat on the notebook to accompany the manuscript:
"Automatic generation of microkinetic mechanisms for heterogeneous catalysis" by
C. Franklin Goldsmith, School of Engineering, Brown University, franklin_goldsmith@brown.edu, and
Richard H. West, Department of Chemical Engineering, Northeastern University, r.west@northeastern.edu

As modified to prepare data for the AIChE conference presentation: 304c Incorporation of Linear Scaling Relations into Automatic Mechanism Generation for Heterogeneous Catalysis Richard H. West, Department of Chemical Engineering, Northeastern University, Boston, MA and C. Franklin Goldsmith, Engineering, Brown University, Providence, RI Tuesday, October 31, 2017 https://aiche.confex.com/aiche/2017/meetingapp.cgi/Paper/500172 https://www.slideshare.net/richardhwest/incorporation-of-linear-scaling-relations-into-automatic-mechanism-generation-for-heterogeneous-catalysis

Then further updated to model catalytic combustion of methanol.

First, we print what git commit we were on when we ran this notebook, for both the source code (RMG-Py) and the database.

In [2]:
%%bash
export RMGpy=~/Code/RMG-Py
pwd
git log -n1 --pretty=oneline
cd ../RMG-database
pwd
git log -n1 --pretty=oneline
/Users/emilymazeau/Code/linear-scaling-tests
d72d26714a6fd8226f70b98333c9e43ae4f9a627 first try to replicate the volcano plot with not regenerated input files
/Users/emilymazeau/Code/RMG-database
48e0a55a7d3f8270c5ffe2dc78f850ed50f7f4f2 Merge remote-tracking branch 'official/master' into cat, again

Model generation¶

We start with a base input file to generate a mechanism for CH3OH plus O2. First we print the input file we'll use to generate the model.

In [3]:
%cat base/input.py
# Data sources
database(
    thermoLibraries=['surfaceThermo', 'primaryThermoLibrary', 'thermo_DFT_CCSDTF12_BAC','DFT_QCI_thermo'],
    reactionLibraries = [('Deutschmann_Ni', False)],
    seedMechanisms = [],
    kineticsDepositories = ['training'],
    kineticsFamilies = 'default',
    kineticsEstimator = 'rate rules',
    bindingEnergies = { # default values for Ni(111)
                       'C':(-5.997, 'eV/molecule'),
                       'H':(-2.778, 'eV/molecule'),
                       'O':(-4.485, 'eV/molecule'),
                       }
)

# List of species
species(
    label='X',
    reactive=True,
    structure=adjacencyList("1 X u0"),
)

species(
    label='CH4',
    reactive=True,
    structure=SMILES("[CH4]"),
)
species(
   label='O2',
   reactive=True,
   structure=adjacencyList(
       """
1 O u1 p2 c0 {2,S}
2 O u1 p2 c0 {1,S}
"""),
)

species(
    label='N2',
    reactive=False,
    structure=SMILES("N#N"),
)

species(
    label='CO2',
    reactive=True,
    structure=SMILES("O=C=O"),
)

species(
    label='H2O',
    reactive=True,
    structure=SMILES("O"),
)

species(
    label='H2',
    reactive=True,
    structure=SMILES("[H][H]"),
)

species(
    label='CO',
    reactive=True,
    structure=SMILES("[C-]#[O+]"),
)

species(
    label='C2H6',
    reactive=True,
    structure=SMILES("CC"),
)

species(
    label='CH2O',
    reactive=True,
    structure=SMILES("C=O"),
)

species(
    label='CH3',
    reactive=True,
    structure=SMILES("[CH3]"),
)

species(
    label='C3H8',
    reactive=True,
    structure=SMILES("CCC"),
)

species(
    label='H',
    reactive=True,
    structure=SMILES("[H]"),
)

species(
    label='C2H5',
    reactive=True,
    structure=SMILES("C[CH2]"),
)

species(
    label='CH3OH',
    reactive=True,
    structure=SMILES("CO"),
)

species(
    label='HCO',
    reactive=True,
    structure=SMILES("[CH]=O"),
)

species(
    label='CH3CHO',
    reactive=True,
    structure=SMILES("CC=O"),
)

species(
    label='OH',
    reactive=True,
    structure=SMILES("[OH]"),
)

species(
    label='C2H4',
    reactive=True,
    structure=SMILES("C=C"),
)

#----------
# Reaction systems
surfaceReactor(
    temperature=(1200,'K'),
    initialPressure=(1.01325, 'bar'),
    initialGasMoleFractions={
        "CH3OH": 0.11764706,
        "O2": 0.17647059,
        "N2": 0.70588235,
    },
    initialSurfaceCoverages={
        "X": 1.0,
    },
    surfaceVolumeRatio=(1.e5, 'm^-1'),
    surfaceSiteDensity=(2.9e-9, 'mol/cm^2'),
#    terminationConversion = { "CH4":0.9,},
    terminationTime=(1.0, 's'),
)

simulator(
    atol=1e-18,
    rtol=1e-12,
)

model(
    toleranceKeepInEdge=0.0,
    toleranceMoveToCore=1e-2,
    toleranceInterruptSimulation=0.1,
    maximumEdgeSpecies=100000
)

options(
    units='si',
    saveRestartPeriod=None,
    generateOutputHTML=True,
    generatePlots=False,
    saveEdgeSpecies=True,
    saveSimulationProfiles=True,
)

Then we try running it. This will take Richard's computer about four minutes and Emily's computer 25 minutes.

In [5]:
%%bash
python ~/Code/RMG-Py/rmg.py base/input.py > /dev/null
tail -n12 base/RMG.log
    Execution time (DD:HH:MM:SS): 00:00:21:34
    Memory used: 1056.74 MB

Making seed mechanism...
Performing final model checks...

MODEL GENERATION COMPLETED

The final model core has 101 species and 324 reactions
The final model edge has 1088 species and 2975 reactions

RMG execution terminated at Thu Jun 21 12:32:54 2018
:root:Removing old /Users/emilymazeau/Code/linear-scaling-tests/base/RMG_backup.log
:root:Moving /Users/emilymazeau/Code/linear-scaling-tests/base/RMG.log to /Users/emilymazeau/Code/linear-scaling-tests/base/RMG_backup.log

/Users/emilymazeau/Code/RMG-Py/rmgpy/data/kinetics/family.py:1135: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.
  data = item.generateReverseRateCoefficient()
/Users/emilymazeau/Code/RMG-Py/rmgpy/thermo/thermoengine.py:58: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.
  wilhoit = thermo0.toWilhoit(B=1000.)
/Users/emilymazeau/.local/lib/python2.7/site-packages/scipy/optimize/optimize.py:1742: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.
  fx = func(x, *args)
/Users/emilymazeau/.local/lib/python2.7/site-packages/scipy/optimize/optimize.py:1794: FutureWarning: `rcond` parameter will change to the default of machine precision times ``max(M, N)`` where M and N are the input matrix dimensions.
To use the future default and silence this warning we advise to pass `rcond=None`, to keep using the old, explicitly pass `rcond=-1`.
  fu = func(x, *args)
/Users/emilymazeau/.local/lib/python2.7/site-packages/matplotlib/cbook/deprecation.py:107: MatplotlibDeprecationWarning: Passing one of 'on', 'true', 'off', 'false' as a boolean is deprecated; use an actual boolean (True/False) instead.
  warnings.warn(message, mplDeprecation, stacklevel=1)

There are 70 species and 218 reactions (?)

Data processing¶

Next we will import some libraries and set things up to start importing and analyzing the simulation results.

In [4]:
%config InlineBackend.figure_format = 'retina'
In [5]:
%matplotlib inline
from matplotlib import pyplot as plt
import matplotlib

# The default output Type 3 (Type3) fonts can't be edited in Adobe Illustrator
# but Type 42 (TrueType) fonts can be, making it easier to move labels around
# slightly to improve layout before publication.
matplotlib.rcParams['pdf.fonttype'] = 42 

# Seaborn helps make matplotlib graphics nicer
import seaborn as sns
sns.set_style("ticks")
sns.set_context("paper", font_scale=1.5, rc={"lines.linewidth": 2.0})

import os
import re
import pandas as pd
import numpy as np
import shutil
import subprocess
import multiprocessing
In [6]:
def get_last_csv_file(job_directory):
    """
    Find the CSV file from the largest model in the provided job directory.
    
    For CSV files named `simulation_1_13.csv` you want 13 to be the highest number.
    """
    solver_directory = os.path.join(job_directory,'solver')
    csv_files = sorted([f for f in os.listdir(solver_directory) if f.endswith('.csv') ],
                       key=lambda f: int(f[:-4].split('_')[2]))
    return os.path.join(solver_directory, csv_files[-1])
    
job_directory = 'base'
get_last_csv_file(job_directory)
Out[6]:
'base/solver/simulation_1_101.csv'

We will use Pandas to import the csv file

In [7]:
last_csv_file = get_last_csv_file(job_directory)
data = pd.read_csv(last_csv_file)
data
Out[7]:
Time (s) Volume (m^3) N2 Ar He Ne X(1) CH4(2) O2(3) CO2(4) ... SX(943) SX(974) C2H3(50) SX(617) C3H8(43) SX(254) C3H7X(41) C3H7X(42) SX(169) C3H7(62)
0 0.000000e+00 0.098469 0.705882 0.0 0.0 0.0 0.285560 0.000000e+00 1.764706e-01 0.000000e+00 ... 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
1 1.856202e-23 0.098469 0.705882 0.0 0.0 0.0 0.285560 4.717598e-50 1.764706e-01 1.680474e-101 ... 1.815842e-157 -9.103578e-194 2.978876e-81 5.589589e-89 -8.244069e-106 -1.239111e-158 -7.465910e-108 -2.488637e-108 -1.757593e-143 -1.956214e-116
2 5.197367e-23 0.098469 0.705882 0.0 0.0 0.0 0.285560 7.129194e-49 1.764706e-01 1.856996e-99 ... 5.575303e-155 1.005338e-190 6.663628e-80 2.893332e-87 1.383778e-104 3.173318e-156 1.253163e-106 4.177209e-107 3.571228e-141 3.283537e-115
3 1.187970e-22 0.098469 0.705882 0.0 0.0 0.0 0.285560 7.669248e-48 1.764706e-01 1.680071e-97 ... 5.325582e-152 4.394961e-187 1.316621e-78 1.263342e-85 1.631208e-102 3.592568e-153 1.477238e-104 4.924127e-105 1.980076e-138 3.870682e-113
4 2.524435e-22 0.098469 0.705882 0.0 0.0 0.0 0.285560 7.049355e-47 1.764706e-01 1.271504e-95 ... 3.638884e-149 1.276425e-183 2.342857e-77 4.703106e-84 1.289394e-100 2.704337e-150 1.167688e-102 3.892294e-103 7.259774e-136 3.059627e-111
5 5.197367e-22 0.098469 0.705882 0.0 0.0 0.0 0.285560 6.030870e-46 1.764706e-01 8.833452e-94 ... 2.142735e-146 3.095097e-180 3.951573e-76 1.620291e-82 9.125113e-99 1.667862e-147 8.263791e-101 2.754597e-101 2.208166e-133 2.165361e-109
6 1.054323e-21 0.098469 0.705882 0.0 0.0 0.0 0.285560 4.986439e-45 1.764706e-01 5.887460e-92 ... 1.175272e-143 6.887802e-177 6.490479e-75 5.377156e-81 6.135928e-97 9.355449e-145 5.556756e-99 1.852252e-99 6.150769e-131 1.456100e-107
7 2.123496e-21 0.098469 0.705882 0.0 0.0 0.0 0.285560 4.054889e-44 1.764706e-01 3.844737e-90 ... 6.226505e-141 1.469950e-173 1.052137e-73 1.752068e-79 4.024417e-95 5.011346e-142 3.644551e-97 1.214850e-97 1.641749e-128 9.551058e-106
8 4.261841e-21 0.098469 0.705882 0.0 0.0 0.0 0.285560 3.270415e-43 1.764706e-01 2.485498e-88 ... 3.242679e-138 3.072856e-170 1.694442e-72 5.657321e-78 2.607249e-93 2.624126e-139 2.361150e-95 7.870500e-96 4.291093e-126 6.188796e-104
9 8.538531e-21 0.098469 0.705882 0.0 0.0 0.0 0.285560 2.626975e-42 1.764706e-01 1.598725e-86 ... 1.674409e-135 6.357963e-167 2.719981e-71 1.818495e-76 1.678830e-91 1.358697e-136 1.520365e-93 5.067884e-94 1.109958e-123 3.986399e-102
10 1.709191e-20 0.098469 0.705882 0.0 0.0 0.0 0.285560 2.105846e-41 1.764706e-01 1.025754e-84 ... 8.609417e-133 1.308786e-163 4.359146e-70 5.832261e-75 1.077724e-89 6.995577e-134 9.759975e-92 3.253325e-92 2.856230e-121 2.560837e-100
11 3.419867e-20 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.686384e-40 1.764706e-01 6.573068e-83 ... 4.417375e-130 2.687252e-160 6.980580e-69 1.868418e-73 6.907926e-88 3.591769e-131 6.255885e-90 2.085295e-90 7.330871e-119 1.643698e-98
12 6.841220e-20 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.349791e-39 1.764706e-01 4.209408e-81 ... 2.264094e-127 5.510526e-157 1.117433e-67 5.982292e-72 4.424433e-86 1.841561e-128 4.006809e-88 1.335603e-88 1.879129e-116 1.055674e-96
13 1.368392e-19 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.080106e-38 1.764706e-01 2.694878e-79 ... 1.159831e-124 1.129277e-153 1.788530e-66 1.914870e-70 2.832713e-84 9.435417e-126 2.565332e-86 8.551107e-87 4.813678e-114 6.796104e-95
14 2.736933e-19 0.098469 0.705882 0.0 0.0 0.0 0.285560 8.641944e-38 1.764706e-01 1.725013e-77 ... 5.939906e-122 2.313497e-150 2.862815e-65 6.128444e-69 1.813280e-82 4.832665e-123 1.642124e-84 5.473748e-85 1.232699e-111 4.397993e-93
15 5.474015e-19 0.098469 0.705882 0.0 0.0 0.0 0.285560 6.913992e-37 1.764706e-01 1.104137e-75 ... 3.041635e-119 4.738796e-147 4.583539e-64 1.961239e-67 1.160610e-80 2.474824e-120 1.051059e-82 3.503531e-83 3.156219e-109 2.875993e-91
16 1.094818e-18 0.098469 0.705882 0.0 0.0 0.0 0.285560 5.531368e-36 1.764706e-01 7.067676e-74 ... 1.557420e-116 9.705825e-144 7.342818e-63 6.276186e-66 7.428253e-79 1.267334e-117 6.727098e-81 2.242366e-81 8.080571e-107 1.918820e-89
17 2.189651e-18 0.098469 0.705882 0.0 0.0 0.0 0.285560 4.425164e-35 1.764706e-01 4.525567e-72 ... 7.974261e-114 1.987830e-140 1.177736e-61 2.008415e-64 4.754193e-77 6.491131e-115 4.305444e-79 1.435148e-79 2.068708e-104 1.328041e-87
18 4.379316e-18 0.098469 0.705882 0.0 0.0 0.0 0.285560 3.540157e-34 1.764706e-01 2.901684e-70 ... 4.082893e-111 4.071150e-137 1.893572e-60 6.426984e-63 3.042718e-75 3.327583e-112 2.755515e-77 9.185051e-78 5.295991e-102 9.779147e-86
19 8.758647e-18 0.098469 0.705882 0.0 0.0 0.0 0.285560 2.832135e-33 1.764706e-01 1.870304e-68 ... 2.090462e-108 8.337765e-134 3.059113e-59 2.056644e-61 1.947349e-73 1.711788e-109 1.763539e-75 5.878462e-76 1.355783e-99 7.896569e-84
20 1.751731e-17 0.098469 0.705882 0.0 0.0 0.0 0.285560 2.265708e-32 1.764706e-01 1.230383e-66 ... 1.070325e-105 1.707569e-130 4.988628e-58 6.581274e-60 1.246305e-71 2.758656e-106 1.128666e-73 3.762219e-74 7.502154e-97 7.150303e-82
21 3.503463e-17 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.812562e-31 1.764706e-01 8.721977e-65 ... 5.480112e-103 3.497053e-127 8.282737e-57 2.106010e-58 7.976329e-70 1.565004e-103 7.223442e-72 2.407814e-72 1.924569e-94 7.259688e-80
22 7.006928e-17 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.450040e-30 1.764706e-01 7.735366e-63 ... 2.805860e-100 7.161743e-124 1.421534e-55 6.739236e-57 5.104812e-68 1.107687e-100 4.622969e-70 1.540989e-70 4.926853e-92 8.081027e-78
23 1.051039e-16 0.098469 0.705882 0.0 0.0 0.0 0.285560 4.304789e-30 1.764706e-01 4.904492e-62 ... 2.412932e-99 7.587052e-123 5.421762e-55 3.143839e-56 2.914548e-67 6.051405e-100 2.639443e-69 8.798144e-70 2.678597e-91 5.270121e-77
24 1.751732e-16 0.098469 0.705882 0.0 0.0 0.0 0.285560 2.016426e-29 1.764706e-01 2.662649e-60 ... 4.863211e-97 6.210470e-120 4.684244e-54 4.212949e-55 7.541711e-66 4.802546e-97 6.829850e-68 2.276616e-68 4.070377e-89 2.238711e-75
25 2.382356e-16 0.098469 0.705882 0.0 0.0 0.0 0.285560 4.653069e-29 1.764706e-01 1.413330e-59 ... 3.609156e-96 6.148633e-119 1.358069e-53 1.398248e-54 3.129736e-65 5.560508e-96 2.834321e-67 9.447736e-68 2.904789e-88 1.025904e-74
26 3.012980e-16 0.098469 0.705882 0.0 0.0 0.0 0.285560 8.959133e-29 1.764706e-01 5.394133e-59 ... 1.516579e-95 3.060522e-118 3.232480e-53 3.757111e-54 9.856511e-65 2.028853e-95 8.926157e-67 2.975385e-67 9.721032e-88 3.701024e-74
27 3.643603e-16 0.098469 0.705882 0.0 0.0 0.0 0.285560 1.523261e-28 1.764706e-01 1.655002e-58 ... 5.483267e-95 1.323297e-117 6.524392e-53 8.311626e-54 2.518899e-64 1.072338e-94 2.281141e-66 7.603801e-67 3.503235e-87 1.042595e-73
28 4.904851e-16 0.098469 0.705882 0.0 0.0 0.0 0.285560 3.367896e-28 1.764706e-01 9.821801e-58 ... 5.328879e-94 2.029880e-116 1.820282e-52 2.615099e-53 1.026388e-63 2.191740e-93 9.295076e-66 3.098358e-66 3.159397e-86 4.781406e-73
29 6.166098e-16 0.098469 0.705882 0.0 0.0 0.0 0.285560 6.309620e-28 1.764706e-01 3.769097e-57 ... 2.428248e-93 1.166697e-115 4.120776e-52 6.468490e-53 3.077094e-63 1.239401e-92 2.786648e-65 9.288826e-66 1.327957e-85 1.575102e-72
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
17892 9.723079e-01 0.098469 0.705882 0.0 0.0 0.0 0.202679 2.440107e-10 2.223801e-12 8.482918e-02 ... 8.841123e-08 3.718231e-17 3.338055e-06 5.720043e-04 2.448712e-11 3.055749e-11 6.341287e-14 1.316198e-13 8.743226e-17 1.219564e-10
17893 9.737086e-01 0.098469 0.705882 0.0 0.0 0.0 0.202693 2.430040e-10 2.222613e-12 8.483527e-02 ... 8.838944e-08 3.716322e-17 3.338941e-06 5.720034e-04 2.448890e-11 3.053279e-11 6.340779e-14 1.316092e-13 8.733135e-17 1.219377e-10
17894 9.751094e-01 0.098469 0.705882 0.0 0.0 0.0 0.202708 2.420029e-10 2.221428e-12 8.484133e-02 ... 8.836771e-08 3.714418e-17 3.339827e-06 5.720025e-04 2.449067e-11 3.050811e-11 6.340272e-14 1.315987e-13 8.723087e-17 1.219191e-10
17895 9.765101e-01 0.098469 0.705882 0.0 0.0 0.0 0.202723 2.410072e-10 2.220246e-12 8.484739e-02 ... 8.834603e-08 3.712519e-17 3.340710e-06 5.720016e-04 2.449245e-11 3.048346e-11 6.339765e-14 1.315882e-13 8.713083e-17 1.219005e-10
17896 9.779108e-01 0.098469 0.705882 0.0 0.0 0.0 0.202738 2.400169e-10 2.219068e-12 8.485343e-02 ... 8.832440e-08 3.710624e-17 3.341592e-06 5.720008e-04 2.449422e-11 3.045882e-11 6.339260e-14 1.315777e-13 8.703122e-17 1.218820e-10
17897 9.790241e-01 0.098469 0.705882 0.0 0.0 0.0 0.202749 2.392337e-10 2.218133e-12 8.485822e-02 ... 8.830725e-08 3.709122e-17 3.342292e-06 5.720001e-04 2.449562e-11 3.043926e-11 6.338860e-14 1.315694e-13 8.695236e-17 1.218673e-10
17898 9.801374e-01 0.098469 0.705882 0.0 0.0 0.0 0.202761 2.384539e-10 2.217200e-12 8.486301e-02 ... 8.829013e-08 3.707624e-17 3.342991e-06 5.719994e-04 2.449702e-11 3.041970e-11 6.338459e-14 1.315611e-13 8.687376e-17 1.218526e-10
17899 9.812506e-01 0.098469 0.705882 0.0 0.0 0.0 0.202772 2.376775e-10 2.216269e-12 8.486779e-02 ... 8.827304e-08 3.706128e-17 3.343689e-06 5.719987e-04 2.449842e-11 3.040016e-11 6.338060e-14 1.315528e-13 8.679544e-17 1.218379e-10
17900 9.823639e-01 0.098469 0.705882 0.0 0.0 0.0 0.202784 2.369045e-10 2.215340e-12 8.487256e-02 ... 8.825599e-08 3.704635e-17 3.344386e-06 5.719980e-04 2.449982e-11 3.038063e-11 6.337661e-14 1.315445e-13 8.671738e-17 1.218233e-10
17901 9.833658e-01 0.098469 0.705882 0.0 0.0 0.0 0.202794 2.362116e-10 2.214506e-12 8.487684e-02 ... 8.824066e-08 3.703294e-17 3.345012e-06 5.719973e-04 2.450108e-11 3.036307e-11 6.337302e-14 1.315371e-13 8.664736e-17 1.218101e-10
17902 9.843678e-01 0.098469 0.705882 0.0 0.0 0.0 0.202805 2.355214e-10 2.213673e-12 8.488112e-02 ... 8.822537e-08 3.701956e-17 3.345638e-06 5.719967e-04 2.450233e-11 3.034551e-11 6.336944e-14 1.315296e-13 8.657755e-17 1.217970e-10
17903 9.853697e-01 0.098469 0.705882 0.0 0.0 0.0 0.202815 2.348339e-10 2.212842e-12 8.488539e-02 ... 8.821010e-08 3.700620e-17 3.346263e-06 5.719961e-04 2.450358e-11 3.032797e-11 6.336587e-14 1.315222e-13 8.650795e-17 1.217839e-10
17904 9.862715e-01 0.098469 0.705882 0.0 0.0 0.0 0.202825 2.342174e-10 2.212095e-12 8.488924e-02 ... 8.819638e-08 3.699420e-17 3.346824e-06 5.719955e-04 2.450471e-11 3.031218e-11 6.336266e-14 1.315155e-13 8.644550e-17 1.217721e-10
17905 9.871732e-01 0.098469 0.705882 0.0 0.0 0.0 0.202834 2.336031e-10 2.211349e-12 8.489307e-02 ... 8.818268e-08 3.698222e-17 3.347385e-06 5.719950e-04 2.450583e-11 3.029641e-11 6.335945e-14 1.315089e-13 8.638321e-17 1.217604e-10
17906 9.880749e-01 0.098469 0.705882 0.0 0.0 0.0 0.202843 2.329909e-10 2.210604e-12 8.489690e-02 ... 8.816900e-08 3.697026e-17 3.347946e-06 5.719944e-04 2.450695e-11 3.028064e-11 6.335624e-14 1.315022e-13 8.632110e-17 1.217486e-10
17907 9.889767e-01 0.098469 0.705882 0.0 0.0 0.0 0.202853 2.323809e-10 2.209861e-12 8.490073e-02 ... 8.815534e-08 3.695831e-17 3.348506e-06 5.719938e-04 2.450808e-11 3.026488e-11 6.335304e-14 1.314956e-13 8.625916e-17 1.217369e-10
17908 9.898784e-01 0.098469 0.705882 0.0 0.0 0.0 0.202862 2.317730e-10 2.209119e-12 8.490455e-02 ... 8.814170e-08 3.694639e-17 3.349065e-06 5.719933e-04 2.450920e-11 3.024913e-11 6.334984e-14 1.314889e-13 8.619739e-17 1.217251e-10
17909 9.907802e-01 0.098469 0.705882 0.0 0.0 0.0 0.202871 2.311673e-10 2.208378e-12 8.490836e-02 ... 8.812808e-08 3.693449e-17 3.349623e-06 5.719927e-04 2.451031e-11 3.023339e-11 6.334665e-14 1.314823e-13 8.613579e-17 1.217134e-10
17910 9.916819e-01 0.098469 0.705882 0.0 0.0 0.0 0.202880 2.305637e-10 2.207638e-12 8.491217e-02 ... 8.811448e-08 3.692260e-17 3.350181e-06 5.719922e-04 2.451143e-11 3.021766e-11 6.334346e-14 1.314757e-13 8.607436e-17 1.217018e-10
17911 9.925837e-01 0.098469 0.705882 0.0 0.0 0.0 0.202890 2.299622e-10 2.206899e-12 8.491598e-02 ... 8.810091e-08 3.691074e-17 3.350739e-06 5.719916e-04 2.451255e-11 3.020193e-11 6.334027e-14 1.314691e-13 8.601309e-17 1.216901e-10
17912 9.934854e-01 0.098469 0.705882 0.0 0.0 0.0 0.202899 2.293627e-10 2.206162e-12 8.491978e-02 ... 8.808735e-08 3.689889e-17 3.351295e-06 5.719911e-04 2.451366e-11 3.018621e-11 6.333709e-14 1.314625e-13 8.595199e-17 1.216784e-10
17913 9.943872e-01 0.098469 0.705882 0.0 0.0 0.0 0.202908 2.287654e-10 2.205426e-12 8.492357e-02 ... 8.807381e-08 3.688706e-17 3.351851e-06 5.719905e-04 2.451478e-11 3.017050e-11 6.333391e-14 1.314559e-13 8.589106e-17 1.216668e-10
17914 9.952889e-01 0.098469 0.705882 0.0 0.0 0.0 0.202917 2.281702e-10 2.204691e-12 8.492736e-02 ... 8.806030e-08 3.687526e-17 3.352407e-06 5.719899e-04 2.451589e-11 3.015480e-11 6.333074e-14 1.314493e-13 8.583030e-17 1.216552e-10
17915 9.961907e-01 0.098469 0.705882 0.0 0.0 0.0 0.202927 2.275771e-10 2.203957e-12 8.493115e-02 ... 8.804680e-08 3.686347e-17 3.352962e-06 5.719894e-04 2.451700e-11 3.013911e-11 6.332757e-14 1.314427e-13 8.576970e-17 1.216435e-10
17916 9.970924e-01 0.098469 0.705882 0.0 0.0 0.0 0.202936 2.269860e-10 2.203224e-12 8.493493e-02 ... 8.803333e-08 3.685170e-17 3.353516e-06 5.719888e-04 2.451811e-11 3.012342e-11 6.332440e-14 1.314361e-13 8.570927e-17 1.216319e-10
17917 9.979942e-01 0.098469 0.705882 0.0 0.0 0.0 0.202945 2.263970e-10 2.202493e-12 8.493870e-02 ... 8.801987e-08 3.683995e-17 3.354070e-06 5.719883e-04 2.451921e-11 3.010775e-11 6.332124e-14 1.314296e-13 8.564900e-17 1.216204e-10
17918 9.988959e-01 0.098469 0.705882 0.0 0.0 0.0 0.202954 2.258100e-10 2.201762e-12 8.494247e-02 ... 8.800643e-08 3.682822e-17 3.354623e-06 5.719877e-04 2.452032e-11 3.009208e-11 6.331808e-14 1.314230e-13 8.558890e-17 1.216088e-10
17919 9.997976e-01 0.098469 0.705882 0.0 0.0 0.0 0.202963 2.252251e-10 2.201033e-12 8.494624e-02 ... 8.799302e-08 3.681651e-17 3.355175e-06 5.719872e-04 2.452143e-11 3.007642e-11 6.331493e-14 1.314165e-13 8.552895e-17 1.215972e-10
17920 1.000000e+00 0.098469 0.705882 0.0 0.0 0.0 0.202965 2.250942e-10 2.200870e-12 8.494708e-02 ... 8.799001e-08 3.681388e-17 3.355299e-06 5.719871e-04 2.452167e-11 3.007290e-11 6.331422e-14 1.314150e-13 8.551553e-17 1.215946e-10
17921 1.000699e+00 0.098469 0.705882 0.0 0.0 0.0 0.202972 2.246423e-10 2.200305e-12 8.495000e-02 ... 8.797962e-08 3.680481e-17 3.355727e-06 5.719866e-04 2.452253e-11 3.006076e-11 6.331177e-14 1.314099e-13 8.546918e-17 1.215857e-10

17922 rows × 103 columns

In [8]:
def get_pandas_data(job_directory):
    """
    Get the last CSV file from the provided job directory,
    import it into a Pandas data frame, and tidy it up a bit.
    """
    last_csv_file = get_last_csv_file(job_directory)
    data = pd.read_csv(last_csv_file)
    
    # Make the Time into an index and remove it as a column
    data.index = data['Time (s)']
    del data['Time (s)']
    # Remove inerts that RMG added automatically but we're not using
    for i in 'Ar He Ne'.split():
        del data[i]
    # Remove the Volume column
    del data['Volume (m^3)']
    # Set any amounts that are slightly negative (within the ODE solver's ATOL tolerance), equal to zero
    # to allow 'area' plots without warnings or errors.
    # Anything more negative than -1e-16 probably indicates a bug in RMG and should not be hidden here ;-)
    data[(data<0) & (data>-1e-16)] = 0
    return data
In [44]:
def rename_columns(data):
    """
    Removes the number (##) from the end of the column names, in place,
    unless doing so would make the names collide.
    Also renames a few species so the plot labels match the names in the manuscript.
    """
    import re
    old = data.columns
    new = [re.sub('\(\d+\)$','',n) for n in old]
    # don't translate names that would no longer be unique
    mapping = {k:v for k,v in zip(old,new) if new.count(v)==1}
    data.rename(columns=mapping, inplace=True)
    
    # Now change a few species that are named differently in the manuscript
    # than in the thermodynamics database used to build the model,
    # so that the plot labels match the manuscript.
    mapping = {'COX':'COvdwX', 'OCX': 'CO=X', 'C2H3X':'CH3CX', 'C2H3OX':'CH3CXO',
               'CO2(4)':'CO2'}
    data.rename(columns=mapping, inplace=True)
In [46]:
data1 = get_pandas_data('base')
rename_columns(data1)
data1.columns
Out[46]:
Index([u'N2', u'X', u'CH4(2)', u'O2', u'CO2', u'H2O(5)', u'H2(6)', u'CO(7)',
       u'C2H6(8)', u'CH2O(9)', u'CH3', u'C3H8(11)', u'H', u'C2H5',
       u'CH3OH(14)', u'HCO', u'CH3CHO(16)', u'OH', u'C2H4(18)', u'CH3OH(46)',
       u'OX', u'HX', u'CH3OX(45)', u'OCCO', u'HOX', u'CH2OX', u'CHOX',
       u'H2(36)', u'H2O(35)', u'HO2(136)', u'SX(156)', u'SX(155)', u'CX',
       u'CO=X', u'CO(37)', u'CH3OX(44)', u'CH3OO', u'CH3X', u'CH2X', u'CHX',
       u'CXHO', u'SX(165)', u'CH4(33)', u'CO2(34)', u'SX(185)', u'SX(183)',
       u'SX(246)', u'SX(225)', u'SX(143)', u'CH2OH', u'HO2(52)', u'CH2O(40)',
       u'HOCXO', u'SX(265)', u'CH3O', u'SX(191)', u'SX(333)', u'SX(356)',
       u'SX(178)', u'HOOH(153)', u'HOOH(408)', u'CH2(S)', u'C2H6O', u'SX(422)',
       u'C2H4(51)', u'C2H4(32)', u'CHO4', u'O(T)', u'CH3O3', u'DME(202)',
       u'DME(494)', u'C3H6O', u'S(439)', u'C2H4O(412)', u'C2H4O(30)',
       u'CH3CHO(49)', u'C2H4O(139)', u'SX(620)', u'S(411)', u'S(692)',
       u'S(693)', u'C2H6(39)', u'CH4O2', u'SX(878)', u'SX(171)', u'SX(885)',
       u'SX(914)', u'SX(161)', u'SX(943)', u'SX(974)', u'C2H3', u'SX(617)',
       u'C3H8(43)', u'SX(254)', u'C3H7X(41)', u'C3H7X(42)', u'SX(169)',
       u'C3H7'],
      dtype='object')
In [47]:
# Test it with some plots
data1[['CH3OH(14)', 'O2']].plot.line()
data1[['CO(7)', 'H2(6)']].plot.line()
data1[['H2O(5)']].plot.line()
Out[47]:
<matplotlib.axes._subplots.AxesSubplot at 0x1165dc390>
In [48]:
species_names = data1.columns
# just the gas phase species that aren't always zero:
gas_phase = [n for n in species_names if 'X' not in n and (data1[n]>0).any()]
# all the surface species
surface_phase = [n for n in species_names if 'X' in n]
surface_phase.remove('X')
data1[gas_phase].plot.line()
data1[surface_phase].plot.line()
Out[48]:
<matplotlib.axes._subplots.AxesSubplot at 0x116cd6710>
In [49]:
print "Significant species (those that exceed 0.001 mol at some point)"
significant = [n for n in data1.columns if(data1[n]>0.001).any()]
with sns.color_palette("hls", len(significant)):
    data1[significant].plot.area(legend='reverse')
Significant species (those that exceed 0.001 mol at some point)
In [50]:
lim = 1e-4
surface = [n for n in data1.columns if 'X' in n and n!='X' and (data1[n]>lim).any() ]
print "The {} surface species that exceed {:.1e} mol at some point".format(len(surface), lim)
total_sites = max(data1['X'])
with sns.color_palette('Set3',len(surface)):
    (data1[surface]/total_sites).plot.area(legend='reverse')
    plt.ylabel('site fraction')
    plt.xlim(0,1e-2) ## notice this is much less than 1 seconde
The 12 surface species that exceed 1.0e-04 mol at some point

Effect of binding energies¶

Now we will use that template 'base' input file to create a ton of other input files with different binding energies, then run them all in RMG-Cat, then process the results.

In [51]:
# First, make a series of input files in separate directories

with open(os.path.join('base', 'input.py')) as infile:
    input_file = infile.read()
    
base_directory = 'binding_energies'
def directory(carbon,oxygen):
    return os.path.join(base_directory, "c{:.3f}o{:.3f}".format(carbon,oxygen))

def make_input(binding_energies):
    """
    Make an input file for the given (carbon,oxygen) tuple (or iterable) of binding energies
    and return the name of the directory in which it is saved.
    """
    carbon, oxygen = binding_energies
    output = input_file
    out_dir = directory(carbon, oxygen)
    carbon_string = "'C':({:f}, 'eV/molecule')".format(carbon)
    output = re.sub("'C':\(.*?, 'eV/molecule'\)", carbon_string, output)
    oxygen_string = "'O':({:f}, 'eV/molecule')".format(oxygen)
    output = re.sub("'O':\(.*?, 'eV/molecule'\)", oxygen_string, output)
    os.path.exists(out_dir) or os.makedirs(out_dir)
    out_file = os.path.join(out_dir, 'input.py')
    with open(out_file,'w') as outfile:
        outfile.write(output)
    shutil.copy(os.path.join('base','run.sh'), out_dir)
    return out_dir

    
print make_input((-8,-3.5))
    
binding_energies/c-8.000o-3.500
In [52]:
def run_simulation(binding_energies):
    """
    Assuming a job file already exists, run it.  This one is local.
    Takes a tuple of binding energies, (carbon, oxygen)
    """
    carbon, oxygen = binding_energies
    job_directory = directory(carbon, oxygen)
    print job_directory
    assert os.path.exists(job_directory)
    return subprocess.check_call('./run.sh', cwd=job_directory)
# a test experiment = (-8,-3.5) make_input(experiment) run_simulation(experiment)
In [53]:
# Revised range
carbon_range = (-7.5, -2)
oxygen_range = (-6.5, -1.5)
grid_size = 9
mesh  = np.mgrid[carbon_range[0]:carbon_range[1]:grid_size*1j, 
                 oxygen_range[0]:oxygen_range[1]:grid_size*1j]

with sns.axes_style("whitegrid"):
    plt.axis('square')
    plt.xlim(carbon_range)
    plt.ylim(oxygen_range)
plt.show()
    
mesh
Out[53]:
array([[[-7.5   , -7.5   , -7.5   , -7.5   , -7.5   , -7.5   , -7.5   ,
         -7.5   , -7.5   ],
        [-6.8125, -6.8125, -6.8125, -6.8125, -6.8125, -6.8125, -6.8125,
         -6.8125, -6.8125],
        [-6.125 , -6.125 , -6.125 , -6.125 , -6.125 , -6.125 , -6.125 ,
         -6.125 , -6.125 ],
        [-5.4375, -5.4375, -5.4375, -5.4375, -5.4375, -5.4375, -5.4375,
         -5.4375, -5.4375],
        [-4.75  , -4.75  , -4.75  , -4.75  , -4.75  , -4.75  , -4.75  ,
         -4.75  , -4.75  ],
        [-4.0625, -4.0625, -4.0625, -4.0625, -4.0625, -4.0625, -4.0625,
         -4.0625, -4.0625],
        [-3.375 , -3.375 , -3.375 , -3.375 , -3.375 , -3.375 , -3.375 ,
         -3.375 , -3.375 ],
        [-2.6875, -2.6875, -2.6875, -2.6875, -2.6875, -2.6875, -2.6875,
         -2.6875, -2.6875],
        [-2.    , -2.    , -2.    , -2.    , -2.    , -2.    , -2.    ,
         -2.    , -2.    ]],

       [[-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ],
        [-6.5   , -5.875 , -5.25  , -4.625 , -4.    , -3.375 , -2.75  ,
         -2.125 , -1.5   ]]])
In [54]:
experiments = mesh.reshape((2,-1)).T
experiments
Out[54]:
array([[-7.5   , -6.5   ],
       [-7.5   , -5.875 ],
       [-7.5   , -5.25  ],
       [-7.5   , -4.625 ],
       [-7.5   , -4.    ],
       [-7.5   , -3.375 ],
       [-7.5   , -2.75  ],
       [-7.5   , -2.125 ],
       [-7.5   , -1.5   ],
       [-6.8125, -6.5   ],
       [-6.8125, -5.875 ],
       [-6.8125, -5.25  ],
       [-6.8125, -4.625 ],
       [-6.8125, -4.    ],
       [-6.8125, -3.375 ],
       [-6.8125, -2.75  ],
       [-6.8125, -2.125 ],
       [-6.8125, -1.5   ],
       [-6.125 , -6.5   ],
       [-6.125 , -5.875 ],
       [-6.125 , -5.25  ],
       [-6.125 , -4.625 ],
       [-6.125 , -4.    ],
       [-6.125 , -3.375 ],
       [-6.125 , -2.75  ],
       [-6.125 , -2.125 ],
       [-6.125 , -1.5   ],
       [-5.4375, -6.5   ],
       [-5.4375, -5.875 ],
       [-5.4375, -5.25  ],
       [-5.4375, -4.625 ],
       [-5.4375, -4.    ],
       [-5.4375, -3.375 ],
       [-5.4375, -2.75  ],
       [-5.4375, -2.125 ],
       [-5.4375, -1.5   ],
       [-4.75  , -6.5   ],
       [-4.75  , -5.875 ],
       [-4.75  , -5.25  ],
       [-4.75  , -4.625 ],
       [-4.75  , -4.    ],
       [-4.75  , -3.375 ],
       [-4.75  , -2.75  ],
       [-4.75  , -2.125 ],
       [-4.75  , -1.5   ],
       [-4.0625, -6.5   ],
       [-4.0625, -5.875 ],
       [-4.0625, -5.25  ],
       [-4.0625, -4.625 ],
       [-4.0625, -4.    ],
       [-4.0625, -3.375 ],
       [-4.0625, -2.75  ],
       [-4.0625, -2.125 ],
       [-4.0625, -1.5   ],
       [-3.375 , -6.5   ],
       [-3.375 , -5.875 ],
       [-3.375 , -5.25  ],
       [-3.375 , -4.625 ],
       [-3.375 , -4.    ],
       [-3.375 , -3.375 ],
       [-3.375 , -2.75  ],
       [-3.375 , -2.125 ],
       [-3.375 , -1.5   ],
       [-2.6875, -6.5   ],
       [-2.6875, -5.875 ],
       [-2.6875, -5.25  ],
       [-2.6875, -4.625 ],
       [-2.6875, -4.    ],
       [-2.6875, -3.375 ],
       [-2.6875, -2.75  ],
       [-2.6875, -2.125 ],
       [-2.6875, -1.5   ],
       [-2.    , -6.5   ],
       [-2.    , -5.875 ],
       [-2.    , -5.25  ],
       [-2.    , -4.625 ],
       [-2.    , -4.    ],
       [-2.    , -3.375 ],
       [-2.    , -2.75  ],
       [-2.    , -2.125 ],
       [-2.    , -1.5   ]])
In [55]:
with sns.axes_style("whitegrid"):
    plt.axis('square')
    plt.xlim(carbon_range)
    plt.ylim(oxygen_range)
    plt.plot(*experiments.T, marker='o', linestyle='none')
In [56]:
map(make_input, experiments)
Out[56]:
['binding_energies/c-7.500o-6.500',
 'binding_energies/c-7.500o-5.875',
 'binding_energies/c-7.500o-5.250',
 'binding_energies/c-7.500o-4.625',
 'binding_energies/c-7.500o-4.000',
 'binding_energies/c-7.500o-3.375',
 'binding_energies/c-7.500o-2.750',
 'binding_energies/c-7.500o-2.125',
 'binding_energies/c-7.500o-1.500',
 'binding_energies/c-6.812o-6.500',
 'binding_energies/c-6.812o-5.875',
 'binding_energies/c-6.812o-5.250',
 'binding_energies/c-6.812o-4.625',
 'binding_energies/c-6.812o-4.000',
 'binding_energies/c-6.812o-3.375',
 'binding_energies/c-6.812o-2.750',
 'binding_energies/c-6.812o-2.125',
 'binding_energies/c-6.812o-1.500',
 'binding_energies/c-6.125o-6.500',
 'binding_energies/c-6.125o-5.875',
 'binding_energies/c-6.125o-5.250',
 'binding_energies/c-6.125o-4.625',
 'binding_energies/c-6.125o-4.000',
 'binding_energies/c-6.125o-3.375',
 'binding_energies/c-6.125o-2.750',
 'binding_energies/c-6.125o-2.125',
 'binding_energies/c-6.125o-1.500',
 'binding_energies/c-5.438o-6.500',
 'binding_energies/c-5.438o-5.875',
 'binding_energies/c-5.438o-5.250',
 'binding_energies/c-5.438o-4.625',
 'binding_energies/c-5.438o-4.000',
 'binding_energies/c-5.438o-3.375',
 'binding_energies/c-5.438o-2.750',
 'binding_energies/c-5.438o-2.125',
 'binding_energies/c-5.438o-1.500',
 'binding_energies/c-4.750o-6.500',
 'binding_energies/c-4.750o-5.875',
 'binding_energies/c-4.750o-5.250',
 'binding_energies/c-4.750o-4.625',
 'binding_energies/c-4.750o-4.000',
 'binding_energies/c-4.750o-3.375',
 'binding_energies/c-4.750o-2.750',
 'binding_energies/c-4.750o-2.125',
 'binding_energies/c-4.750o-1.500',
 'binding_energies/c-4.062o-6.500',
 'binding_energies/c-4.062o-5.875',
 'binding_energies/c-4.062o-5.250',
 'binding_energies/c-4.062o-4.625',
 'binding_energies/c-4.062o-4.000',
 'binding_energies/c-4.062o-3.375',
 'binding_energies/c-4.062o-2.750',
 'binding_energies/c-4.062o-2.125',
 'binding_energies/c-4.062o-1.500',
 'binding_energies/c-3.375o-6.500',
 'binding_energies/c-3.375o-5.875',
 'binding_energies/c-3.375o-5.250',
 'binding_energies/c-3.375o-4.625',
 'binding_energies/c-3.375o-4.000',
 'binding_energies/c-3.375o-3.375',
 'binding_energies/c-3.375o-2.750',
 'binding_energies/c-3.375o-2.125',
 'binding_energies/c-3.375o-1.500',
 'binding_energies/c-2.688o-6.500',
 'binding_energies/c-2.688o-5.875',
 'binding_energies/c-2.688o-5.250',
 'binding_energies/c-2.688o-4.625',
 'binding_energies/c-2.688o-4.000',
 'binding_energies/c-2.688o-3.375',
 'binding_energies/c-2.688o-2.750',
 'binding_energies/c-2.688o-2.125',
 'binding_energies/c-2.688o-1.500',
 'binding_energies/c-2.000o-6.500',
 'binding_energies/c-2.000o-5.875',
 'binding_energies/c-2.000o-5.250',
 'binding_energies/c-2.000o-4.625',
 'binding_energies/c-2.000o-4.000',
 'binding_energies/c-2.000o-3.375',
 'binding_energies/c-2.000o-2.750',
 'binding_energies/c-2.000o-2.125',
 'binding_energies/c-2.000o-1.500']
In [57]:
# Now run the simulations using a pool.
# Don't run this cell unless you have a while to wait!! ###

#pool = multiprocessing.Pool()
#result = pool.map(run_simulation, experiments)

# Instead, upload the input files to a cluster and run `start_all.sh`
# The script `stage_all.sh` will stage (add) the results in git
# so you can commit and get them back to analyze with the subsequent cells:
In [58]:
base_directory = 'binding_energies'
# base_directory = 'binding_energies_local'
In [59]:
def get_data(experiment):
    carbon, oxygen = experiment
    directory(carbon,oxygen)
    try:
        data = get_pandas_data(directory(carbon,oxygen))
    except OSError:
        return None
    rename_columns(data)
    return data
In [61]:
datas = {tuple(e): get_data(e) for e in experiments}
In [63]:
def get_max_co2(experiment):
    data = datas[tuple(experiment)]
    if data is None:
        return np.nan
    return data[['CO2']].max()
highest_co2 = max([float(get_max_co2(e)) for e in experiments])

Fit an exponential¶

For this first attempt to extract "rate" we fit an exponential growth curve to the normalized CO2 concentration profile of each simulation. First we plot all the curves on one plot, then we'll plot each with its fitted exponential.

In [64]:
import seaborn as sns
plt.figure(figsize=(5, 4))
num_lines = len(experiments)
sns.set_palette(sns.color_palette("coolwarm",num_lines))

ax = plt.axes()
def make_label(experiment):
    return "{:+.1f}, {:+.1f}".format(*experiment)

for experiment in experiments:
    #print experiment
    data = get_data(experiment)
    if data is None:
        print("No data for {}".format(experiment))
        continue
    times = np.array(data.index)
    normalized = data[['CO2']].values / highest_co2
    ax.plot(times, normalized, label=make_label(experiment))
    #normalized.plot.line(ax=ax, label=make_label(experiment))
    #ax.text(times[-10], normalized[-10], make_label(experiment))

#plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.xlabel('Time (s)')
plt.ylabel('Normalized CO$_2$ concentration')
Out[64]:
Text(0,0.5,'Normalized CO$_2$ concentration')
In [65]:
sns.set_palette('Set1')

def my_function(time, rate):
    "Thing we want to fit."
    return 1. - np.exp(-1*time*rate)
import scipy.optimize


def fit_rate(data):
    normalized = data[['CO2']] / highest_co2
    x_data = np.array(normalized.index)
    y_data = normalized.values[:,0] 
    
    try:
        popt, pcov = scipy.optimize.curve_fit(my_function, x_data, y_data)
    except RuntimeError as e:
        print("☝️ couldn't fit: {}".format(e.message))
        return np.nan
        
    optimal_parameters = popt
    parameter_errors = np.sqrt(np.diag(pcov))
    print("Rate: {} +/- {} (1 st. dev.)".format(optimal_parameters[0],parameter_errors[0]))

    plt.plot(x_data, y_data, 'o')
    plt.plot(x_data, my_function(x_data, *optimal_parameters))
    plt.xlabel('Time (s)')
    plt.ylabel('Normalized CO$_2$ concentration')
    
    plt.show()
    fitted_rate = optimal_parameters[0]
    return fitted_rate

rates = []
for experiment in experiments:
    print experiment
    data = get_data(experiment)
    if data is None:
        print("☝️No data")
        rates.append(np.nan)
        continue
    rate = fit_rate(data)
    rates.append(rate)

rates
    
[-7.5 -6.5]
Rate: 0.0116283698207 +/- 0.000108932362375 (1 st. dev.)
[-7.5   -5.875]
Rate: 0.269756647562 +/- 0.00104801721179 (1 st. dev.)
[-7.5  -5.25]
Rate: 1.24937115277 +/- 0.025492925756 (1 st. dev.)
[-7.5   -4.625]
Rate: 301.208704472 +/- 1.76992922927 (1 st. dev.)
[-7.5 -4. ]
Rate: 818.906522293 +/- 2.70814372908 (1 st. dev.)
[-7.5   -3.375]
Rate: 204.14904304 +/- 2.13256924428 (1 st. dev.)
[-7.5  -2.75]
Rate: 4.02018320404 +/- 0.0586668871366 (1 st. dev.)
[-7.5   -2.125]
Rate: 1908.96249249 +/- 20.9353317166 (1 st. dev.)
[-7.5 -1.5]
Rate: 477.233556901 +/- 5.26086818625 (1 st. dev.)
[-6.8125 -6.5   ]
Rate: 0.00844562607958 +/- 4.10147634019e-05 (1 st. dev.)
[-6.8125 -5.875 ]
Rate: 0.232735848154 +/- 0.00115727772251 (1 st. dev.)
[-6.8125 -5.25  ]
Rate: 13.7543046179 +/- 0.137908819542 (1 st. dev.)
[-6.8125 -4.625 ]
Rate: 42.7801823968 +/- 0.222514330102 (1 st. dev.)
[-6.8125 -4.    ]
Rate: 677.479690608 +/- 2.29444801108 (1 st. dev.)
[-6.8125 -3.375 ]
Rate: 186.452644353 +/- 1.59000928896 (1 st. dev.)
[-6.8125 -2.75  ]
Rate: 7.22630880306 +/- 0.105336886247 (1 st. dev.)
[-6.8125 -2.125 ]
Rate: 1865.90877671 +/- 22.8257642769 (1 st. dev.)
[-6.8125 -1.5   ]
Rate: 264.323315902 +/- 3.2401646453 (1 st. dev.)
[-6.125 -6.5  ]
Rate: 0.00325074019785 +/- 1.36307744148e-05 (1 st. dev.)
[-6.125 -5.875]
Rate: 0.0933392418477 +/- 0.000389002630569 (1 st. dev.)
[-6.125 -5.25 ]
Rate: 2.56420631795 +/- 0.0141393289206 (1 st. dev.)
[-6.125 -4.625]
Rate: 1.28391097513 +/- 0.00826902488831 (1 st. dev.)
[-6.125 -4.   ]
Rate: 121.42109308 +/- 0.878260995559 (1 st. dev.)
[-6.125 -3.375]
Rate: 137.378119742 +/- 1.11918208434 (1 st. dev.)
[-6.125 -2.75 ]
Rate: 4.35132811849 +/- 0.0459997827184 (1 st. dev.)
[-6.125 -2.125]
Rate: 1.45326214186 +/- 0.0347361857237 (1 st. dev.)
[-6.125 -1.5  ]
Rate: 243.079812761 +/- 3.28876394853 (1 st. dev.)
[-5.4375 -6.5   ]
Rate: 4.65755629004e-05 +/- 7.41809017976e-07 (1 st. dev.)
[-5.4375 -5.875 ]
Rate: 0.000524867301661 +/- 9.98610874236e-06 (1 st. dev.)
[-5.4375 -5.25  ]
Rate: 0.00190760580975 +/- 2.30954098514e-05 (1 st. dev.)
[-5.4375 -4.625 ]
Rate: 0.00162635581426 +/- 1.57693931581e-05 (1 st. dev.)
[-5.4375 -4.    ]
Rate: 4.41117762681 +/- 0.0382081196939 (1 st. dev.)
[-5.4375 -3.375 ]
Rate: 67.4276823555 +/- 0.687121463441 (1 st. dev.)
[-5.4375 -2.75  ]
Rate: 4.9405491284 +/- 0.0432924047441 (1 st. dev.)
[-5.4375 -2.125 ]
Rate: 1.30084203542 +/- 0.0207642227641 (1 st. dev.)
[-5.4375 -1.5   ]
Rate: 140.380915666 +/- 1.69130791666 (1 st. dev.)
[-4.75 -6.5 ]
Rate: 8.85058631795e-06 +/- 1.7571138564e-07 (1 st. dev.)
[-4.75  -5.875]
Rate: 1.88086123906e-05 +/- 3.77360054118e-07 (1 st. dev.)
[-4.75 -5.25]
Rate: 1.83269966703e-05 +/- 3.36677705467e-07 (1 st. dev.)
[-4.75  -4.625]
Rate: 1.76552611629e-05 +/- 2.88771097219e-07 (1 st. dev.)
[-4.75 -4.  ]
Rate: 0.423420668226 +/- 0.00426632884957 (1 st. dev.)
[-4.75  -3.375]
Rate: 1.15432707332 +/- 0.0179224331075 (1 st. dev.)
[-4.75 -2.75]
Rate: 1.06609041838 +/- 0.0236561357363 (1 st. dev.)
[-4.75  -2.125]
Rate: 2.97805625337 +/- 0.0643419702822 (1 st. dev.)
[-4.75 -1.5 ]
Rate: 5.29696107701 +/- 0.0218644626509 (1 st. dev.)
[-4.0625 -6.5   ]
Rate: 4.30421204604e-07 +/- 9.74703605462e-09 (1 st. dev.)
[-4.0625 -5.875 ]
Rate: 5.90785218174e-07 +/- 1.41217069991e-08 (1 st. dev.)
[-4.0625 -5.25  ]
Rate: 5.3261698964e-07 +/- 1.18894275176e-08 (1 st. dev.)
[-4.0625 -4.625 ]
Rate: 3.5823425461e-07 +/- 7.58637128e-09 (1 st. dev.)
[-4.0625 -4.    ]
Rate: 0.225442424212 +/- 0.0011590133072 (1 st. dev.)
[-4.0625 -3.375 ]
Rate: 0.323027624227 +/- 0.00265952103219 (1 st. dev.)
[-4.0625 -2.75  ]
Rate: 0.827838216428 +/- 0.0114314433683 (1 st. dev.)
[-4.0625 -2.125 ]
Rate: 0.288420839102 +/- 0.00337408511387 (1 st. dev.)
[-4.0625 -1.5   ]
Rate: 0.489282388109 +/- 0.00181362947902 (1 st. dev.)
[-3.375 -6.5  ]
Rate: 7.25281366801e-09 +/- 8.71844934333e-11 (1 st. dev.)
[-3.375 -5.875]
/Users/emilymazeau/.local/lib/python2.7/site-packages/scipy/optimize/minpack.py:785: OptimizeWarning: Covariance of the parameters could not be estimated
  category=OptimizeWarning)
Rate: -1.77022938017e-10 +/- inf (1 st. dev.)
[-3.375 -5.25 ]
Rate: -1.36519361872e-10 +/- inf (1 st. dev.)
[-3.375 -4.625]
Rate: -4.57460872004e-11 +/- inf (1 st. dev.)
[-3.375 -4.   ]
Rate: 0.0379593936683 +/- 0.000235749016969 (1 st. dev.)
[-3.375 -3.375]
Rate: 0.546394958974 +/- 0.00210204588082 (1 st. dev.)
[-3.375 -2.75 ]
Rate: 0.604168676525 +/- 0.00351642957827 (1 st. dev.)
[-3.375 -2.125]
Rate: 0.0515997904529 +/- 0.000911379742891 (1 st. dev.)
[-3.375 -1.5  ]
Rate: 0.00758278950855 +/- 3.20596668764e-05 (1 st. dev.)
[-2.6875 -6.5   ]
Rate: -1.06168612797e-10 +/- inf (1 st. dev.)
[-2.6875 -5.875 ]
Rate: -2.20065432087e-10 +/- inf (1 st. dev.)
[-2.6875 -5.25  ]
Rate: -7.12002106406e-10 +/- inf (1 st. dev.)
[-2.6875 -4.625 ]
Rate: -2.59030891367e-10 +/- inf (1 st. dev.)
[-2.6875 -4.    ]
Rate: 0.01349074127 +/- 6.90261347363e-05 (1 st. dev.)
[-2.6875 -3.375 ]
Rate: 0.331853826921 +/- 0.000321066530981 (1 st. dev.)
[-2.6875 -2.75  ]
Rate: 0.13826603044 +/- 0.000642846703333 (1 st. dev.)
[-2.6875 -2.125 ]
Rate: 0.0104565104277 +/- 0.000203727072218 (1 st. dev.)
[-2.6875 -1.5   ]
Rate: 0.00180713063762 +/- 3.0355993113e-06 (1 st. dev.)
[-2.  -6.5]
/Users/emilymazeau/anaconda2/envs/rmg_env/lib/python2.7/site-packages/ipykernel_launcher.py:17: DeprecationWarning: BaseException.message has been deprecated as of Python 2.6
☝️ couldn't fit: Optimal parameters not found: Number of calls to function has reached maxfev = 400.
[-2.    -5.875]
☝️ couldn't fit: Optimal parameters not found: Number of calls to function has reached maxfev = 400.
[-2.   -5.25]
☝️ couldn't fit: Optimal parameters not found: Number of calls to function has reached maxfev = 400.
[-2.    -4.625]
☝️ couldn't fit: Optimal parameters not found: Number of calls to function has reached maxfev = 400.
[-2. -4.]
Rate: 0.00219478995249 +/- 1.90046792488e-05 (1 st. dev.)
[-2.    -3.375]
Rate: 0.0244154728466 +/- 3.41376464246e-06 (1 st. dev.)
[-2.   -2.75]
Rate: 0.0016460292398 +/- 1.04634844181e-05 (1 st. dev.)
[-2.    -2.125]
Rate: -6.75197411702e-11 +/- inf (1 st. dev.)
[-2.  -1.5]
Rate: 0.00136887287413 +/- 5.33005436796e-06 (1 st. dev.)
Out[65]:
[0.011628369820670437,
 0.2697566475617567,
 1.2493711527677058,
 301.2087044720744,
 818.9065222928884,
 204.14904304036077,
 4.020183204042752,
 1908.9624924931074,
 477.2335569011,
 0.008445626079579075,
 0.23273584815398152,
 13.754304617875793,
 42.78018239677485,
 677.4796906084193,
 186.4526443530116,
 7.226308803063021,
 1865.9087767109027,
 264.3233159016064,
 0.003250740197848818,
 0.09333924184771168,
 2.5642063179520913,
 1.2839109751334936,
 121.42109308004345,
 137.37811974230988,
 4.351328118486083,
 1.453262141857439,
 243.07981276052973,
 4.65755629003896e-05,
 0.0005248673016612709,
 0.0019076058097541406,
 0.0016263558142569489,
 4.411177626811544,
 67.42768235545174,
 4.940549128400247,
 1.300842035418964,
 140.38091566578214,
 8.850586317948531e-06,
 1.88086123905906e-05,
 1.8326996670264623e-05,
 1.765526116289123e-05,
 0.4234206682258973,
 1.1543270733205089,
 1.0660904183842688,
 2.9780562533680057,
 5.296961077012625,
 4.304212046041523e-07,
 5.907852181738881e-07,
 5.326169896397277e-07,
 3.582342546103414e-07,
 0.2254424242120894,
 0.323027624226705,
 0.8278382164284473,
 0.288420839102247,
 0.48928238810892605,
 7.252813668006008e-09,
 -1.7702293801660128e-10,
 -1.3651936187222933e-10,
 -4.5746087200444616e-11,
 0.03795939366825307,
 0.5463949589737885,
 0.6041686765249897,
 0.05159979045288473,
 0.007582789508547501,
 -1.0616861279680241e-10,
 -2.2006543208701816e-10,
 -7.120021064063521e-10,
 -2.5903089136668737e-10,
 0.013490741270039025,
 0.3318538269211214,
 0.13826603043993463,
 0.010456510427725005,
 0.0018071306376200005,
 nan,
 nan,
 nan,
 nan,
 0.002194789952490837,
 0.024415472846566845,
 0.0016460292397990647,
 -6.751974117020282e-11,
 0.0013688728741255794]
In [66]:
rates = np.array(rates)
fixed_rates = rates * (rates>0) + (1e-9 * (rates<0))
log_rates = np.log(fixed_rates)
log_rates
/Users/emilymazeau/anaconda2/envs/rmg_env/lib/python2.7/site-packages/ipykernel_launcher.py:2: RuntimeWarning: invalid value encountered in greater
  
/Users/emilymazeau/anaconda2/envs/rmg_env/lib/python2.7/site-packages/ipykernel_launcher.py:2: RuntimeWarning: invalid value encountered in less
  
Out[66]:
array([ -4.45430749,  -1.31023503,   0.22264035,   5.70780339,
         6.70796994,   5.31885033,   1.39132747,   7.55431518,
         6.16800601,  -4.7741066 ,  -1.45785117,   2.62135184,
         3.75607497,   6.51837958,   5.22817729,   1.97772837,
         7.53150349,   5.57717304,  -5.72887256,  -2.37151466,
         0.941649  ,   0.24991087,   4.79926461,   4.92273712,
         1.47048111,   0.37381078,   5.49338984,  -9.97443456,
        -7.55236509,  -6.26190633,  -6.42141346,   1.48414169,
         4.21105565,   1.59747648,   0.26301177,   4.94435955,
       -11.63502685, -10.88119569, -10.90713536, -10.94447674,
        -0.85938911,   0.14351755,   0.06399814,   1.09127082,
         1.66713327, -14.65850156, -14.34181331, -14.44546326,
       -14.84207872,  -1.48969048,  -1.13001744,  -0.18893753,
        -1.24333462,  -0.71481548, -18.74187635, -20.72326584,
       -20.72326584, -20.72326584,  -3.27123828,  -0.6044132 ,
        -0.50390185,  -2.96423767,  -4.88187414, -20.72326584,
       -20.72326584, -20.72326584, -20.72326584,  -4.30575166,
        -1.10306069,  -1.97857569,  -4.56053049,  -6.31601497,
                nan,          nan,          nan,          nan,
        -6.12166893,  -3.71253821,  -6.40938941, -20.72326584,
        -6.5937676 ])
In [67]:
rate_grid = np.reshape(log_rates, (grid_size,grid_size))
In [68]:
ax = sns.heatmap(rate_grid)
In [69]:
extent = carbon_range + oxygen_range
extent
Out[69]:
(-7.5, -2, -6.5, -1.5)
In [70]:
# Because the center of a corner pixel is in fact the corner of the grid
# we want to stretch the image a little
c_step = mesh[0,1,0]-mesh[0,0,0]
o_step = mesh[1,0,1]-mesh[1,0,0]
carbon_range2 = (carbon_range[0]-c_step/2, carbon_range[1]+c_step/2)
oxygen_range2 = (oxygen_range[0]-c_step/2, oxygen_range[1]+c_step/2)
extent2 = carbon_range2 + oxygen_range2
extent2
Out[70]:
(-7.84375, -1.65625, -6.84375, -1.15625)
In [71]:
plt.imshow(rate_grid.T, interpolation='none', origin='lower', extent=extent2, aspect='equal')
plt.plot(-5.997, -4.485, 'ok')
plt.text(-5.997, -4.485, 'Ni(111)')
Out[71]:
Text(-5.997,-4.485,'Ni(111)')
In [72]:
# Binding energies extracted from:
u"""
Medford, A. J.; Lausche, A. C.; Abild-Pedersen, F.;
Temel, B.; Schjødt, N. C.; Nørskov, J. K.; Studt, F. 
Activity and Selectivity Trends in 
Synthesis Gas Conversion to Higher Alcohols. 
Topics in Catalysis 2014, 57 (1-4), 135–142 
DOI: 10.1007/s11244-013-0169-0.
"""
medford_energies = { # Carbon, then Oxygen
'Ru': ( 0.010349288486416697, -2.8153856448231256),
'Rh': ( 0.16558861578266493, -2.546620868091181),
'Ni': ( 0.3001293661060802, -2.5881741535441853),
'Ir': ( 0.36222509702457995, -2.826185484230718),
'Pd': ( 0.28460543337645516, -1.207119596734621),
'Pt': ( 0.8796895213454077, -1.445820136503547),
'Cu': ( 2.323415265200518, -1.7218249542757729),
'Ag': ( 3.855109961190168, -0.8341504215550701),
'Au': ( 3.5601552393272975, -0.10963108355266138),
}
In [73]:
# Shift Medford's energies so that Ni matches Wayne Blaylock's Ni
blaylock_ni = np.array([-5.997, -4.485])
old_ni = np.array(medford_energies['Ni'])
shifted_energies = {metal: tuple(blaylock_ni + np.array(E)-old_ni) for metal,E in medford_energies.items()}
shifted_energies
Out[73]:
{'Ag': (-2.442019404915912, -2.730976268010885),
 'Au': (-2.7369741267787826, -2.006456930008476),
 'Cu': (-3.973714100905562, -3.618650800731588),
 'Ir': (-5.9349042690815, -4.723011330686532),
 'Ni': (-5.997, -4.484999999999999),
 'Pd': (-6.012523932729625, -3.1039454431904363),
 'Pt': (-5.417439844760672, -3.3426459829593616),
 'Rh': (-6.131540750323415, -4.443446714546996),
 'Ru': (-6.286780077619664, -4.712211491278941)}
In [74]:
plt.imshow(rate_grid.T, interpolation='none', origin='lower', extent=extent2, aspect='equal')
for metal, coords in shifted_energies.iteritems():
    plt.plot(coords[0], coords[1], 'ok')
    plt.text(coords[0], coords[1], metal)
plt.xlim(carbon_range)
plt.ylim(oxygen_range)
Out[74]:
(-6.5, -1.5)
In [75]:
# Binding energies for close packed surfaces from
"""
Abild-Pedersen, F.; Greeley, J.; Studt, F.; Rossmeisl, J.;
Munter, T. R.; Moses, P. G.; Skúlason, E.; Bligaard, T.; Norskov, J. K. 
Scaling Properties of Adsorption Energies for Hydrogen-Containing 
Molecules on Transition-Metal Surfaces. 
Phys. Rev. Lett. 2007, 99 (1), 016105 
DOI: 10.1103/PhysRevLett.99.016105.
"""
abildpedersen_energies = { # Carbon, then Oxygen
'Ru': ( -6.397727272727272, -5.104763568600047),
'Rh': ( -6.5681818181818175, -4.609771721406942),
'Ni': ( -6.045454545454545, -4.711681807593758),
'Ir': ( -6.613636363636363, -5.94916142557652),
'Pd': ( -6, -3.517877940833916),
'Pt': ( -6.363636363636363, -3.481481481481482),
'Cu': ( -4.159090909090907, -3.85272536687631),
'Ag': ( -2.9545454545454533, -2.9282552993244817),
'Au': ( -3.7499999999999973, -2.302236198462614),
}
In [76]:
abildpedersen_energies['Pt'][0] - abildpedersen_energies['Ni'][0]
Out[76]:
-0.31818181818181834
In [77]:
plt.imshow(rate_grid.T, interpolation='none', origin='lower', extent=extent2, aspect='equal')
for metal, coords in abildpedersen_energies.iteritems():
    color = {'Ag':'w','Au':'w','Cu':'w'}.get(metal,'k')
    plt.plot(coords[0], coords[1], 'o'+color)
    plt.text(coords[0], coords[1], metal, color=color)
plt.xlim(carbon_range)
plt.ylim(oxygen_range)
plt.xlabel('$\Delta E^C$ (eV)')
plt.ylabel('$\Delta E^O$ (eV)')
Out[77]:
Text(0,0.5,'$\\Delta E^O$ (eV)')

Alternative rate definition¶

For this, we plot the CO2 concentration profiles on a log-log plot, and see the characteristic time as the time at which it crosses the diagonal, i.e. when does $\left( \log_{10}(time) + \log_{10}(\frac{CO_2}{CO_2max}) \right) \ge -10$

In [78]:
"""
Plot everything on a log-log plot
"""
import seaborn as sns
plt.figure(figsize=(5, 4))
num_lines = len(experiments)
sns.set_palette(sns.color_palette("coolwarm",num_lines))
ax = plt.axes()

def make_label(experiment):
    return "{:+.1f}, {:+.1f}".format(*experiment)

for experiment in experiments:
    print experiment
    data = get_data(experiment)
    if data is None:
        print("☝️No data")
        continue
    times = np.array(data.index)
    normalized = data[['CO2']].values[:,0] / highest_co2
    ax.loglog(times, normalized, label=make_label(experiment))
    try:
        i = (np.nonzero((np.log10(times)+np.log10(normalized)) > -10))[0][0]
    except IndexError:
        i = -1
    print i, times[i], normalized[i]
    ax.text(times[i], normalized[i], make_label(experiment), rotation=45, ha='center', va='center', fontsize=12)
    
plt.plot([1e-10, 1],[1,1e-10], 'g:')
plt.ylim(1e-10, 1)
plt.xlim(1e-10, 1)
#plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.xlabel('Time (s)')
plt.ylabel('Normalized CO2 concentration')
[-7.5 -6.5]
/Users/emilymazeau/anaconda2/envs/rmg_env/lib/python2.7/site-packages/ipykernel_launcher.py:23: RuntimeWarning: divide by zero encountered in log10
6045 0.0007209773643256943 1.4043825580947508e-07
[-7.5   -5.875]
4430 2.3631536530199522e-05 4.2674521522443015e-06
[-7.5  -5.25]
3776 1.615449449287374e-06 6.241321163788782e-05
[-7.5   -4.625]
2819 1.2315415655980443e-06 8.123931823063947e-05
[-7.5 -4. ]
2737 9.77301102457786e-07 0.00010437820990665065
[-7.5   -3.375]
2523 2.96826638642026e-07 0.0003424375444024062
[-7.5  -2.75]
2253 2.459645540651213e-07 0.00040743918795815436
[-7.5   -2.125]
2873 2.5760573649343595e-07 0.00038852554261268875
[-7.5 -1.5]
2828 4.809285021475257e-07 0.00020807196706713345
[-6.8125 -6.5   ]
5025 0.00019770939221776848 5.084655146809123e-07
[-6.8125 -5.875 ]
2989 6.717351516407535e-06 1.4975800848216095e-05
[-6.8125 -5.25  ]
3916 1.9054661205962455e-06 5.2766286840678355e-05
[-6.8125 -4.625 ]
3438 2.139646198025301e-06 4.6821140515223193e-05
[-6.8125 -4.    ]
3167 1.1943344177101109e-06 8.439086479340919e-05
[-6.8125 -3.375 ]
2619 3.830500056478119e-07 0.00026378794838804277
[-6.8125 -2.75  ]
2226 3.603743019062733e-07 0.0002787956852393076
[-6.8125 -2.125 ]
2567 3.821333498579415e-07 0.00026522935678660115
[-6.8125 -1.5   ]
2476 1.355168623265055e-06 7.422366954419431e-05
[-6.125 -6.5  ]
5404 0.00012534228258049218 7.991394833481019e-07
[-6.125 -5.875]
2896 7.011316330436784e-06 1.427252128036132e-05
[-6.125 -5.25 ]
3736 5.2676051326400975e-06 1.903202379962261e-05
[-6.125 -4.625]
3448 6.0268929031600825e-06 1.664126942229585e-05
[-6.125 -4.   ]
2875 1.944938198660659e-06 5.141861545672312e-05
[-6.125 -3.375]
2881 6.484230819200632e-07 0.0001555733277211437
[-6.125 -2.75 ]
2396 6.393697240676095e-07 0.00015755671846208705
[-6.125 -2.125]
2561 7.907505854855156e-07 0.00012729222401564275
[-6.125 -1.5  ]
2327 4.772661030995205e-06 2.105741896713845e-05
[-5.4375 -6.5   ]
5086 0.00018029655547746002 5.551228079087138e-07
[-5.4375 -5.875 ]
2934 1.956230953898686e-05 5.125288628461189e-06
[-5.4375 -5.25  ]
3379 2.065573869850726e-05 4.8639503051495115e-06
[-5.4375 -4.625 ]
3588 2.492441964209076e-05 4.013823123817657e-06
[-5.4375 -4.    ]
3305 3.5133067129763e-06 2.849713031576345e-05
[-5.4375 -3.375 ]
2771 1.1341017715415605e-06 8.829479541595529e-05
[-5.4375 -2.75  ]
2720 1.181564107299237e-06 8.545911247736673e-05
[-5.4375 -2.125 ]
2603 2.526172624118653e-06 3.97343000055191e-05
[-5.4375 -1.5   ]
2504 1.1523728020874843e-05 8.678197875197882e-06
[-4.75 -6.5 ]
3813 0.00019725789681378258 5.073238344351607e-07
[-4.75  -5.875]
3376 9.806521591746093e-05 1.0200501428965693e-06
[-4.75 -5.25]
3927 0.0001239811447173297 8.106565818468363e-07
[-4.75  -4.625]
4104 0.00014353380474816483 7.004865662027675e-07
[-4.75 -4.  ]
3292 8.995082730659164e-06 1.1117859645969405e-05
[-4.75  -3.375]
2846 2.737245199863232e-06 3.66693349495933e-05
[-4.75 -2.75]
3102 3.2710978340136706e-06 3.0757072124888876e-05
[-4.75  -2.125]
2854 4.910189712195749e-06 2.0617827116682422e-05
[-4.75 -1.5 ]
2482 2.0838363316624526e-05 4.822949869723321e-06
[-4.0625 -6.5   ]
4819 0.0011558998685879506 8.657276739467334e-08
[-4.0625 -5.875 ]
4333 0.0011185130422842672 8.942862253130947e-08
[-4.0625 -5.25  ]
4707 0.0015398719393465619 6.500021561037467e-08
[-4.0625 -4.625 ]
5199 0.0019217650370103871 5.2256196215368666e-08
[-4.0625 -4.    ]
3642 4.921876253569548e-05 2.0379424336575903e-06
[-4.0625 -3.375 ]
4377 1.211608047086975e-05 8.279702509891634e-06
[-4.0625 -2.75  ]
2445 7.74941212945325e-06 1.321976770188403e-05
[-4.0625 -2.125 ]
3142 7.137676549287648e-06 1.4043191343082034e-05
[-4.0625 -1.5   ]
3165 4.478264224206601e-05 2.2547470671725414e-06
[-3.375 -6.5  ]
9376 0.04728210065213861 2.116287443835175e-09
[-3.375 -5.875]
9358 0.04965392097844989 2.0157122804346173e-09
[-3.375 -5.25 ]
9214 0.0714823469758213 1.3997322825013118e-09
[-3.375 -4.625]
9661 0.09089881092241776 1.1007899676379092e-09
[-3.375 -4.   ]
4072 0.00021925284030170343 4.5772700165646984e-07
[-3.375 -3.375]
3692 4.659446827061543e-05 2.1566828367732276e-06
[-3.375 -2.75 ]
2657 8.971529921114449e-06 1.150033277610266e-05
[-3.375 -2.125]
2748 1.074845485267968e-05 9.634106703321613e-06
[-3.375 -1.5  ]
3137 8.522366460079624e-05 1.1737483456038348e-06
[-2.6875 -6.5   ]
-1 1.0003043893782788 8.129462581120606e-11
[-2.6875 -5.875 ]
-1 1.0022276349201975 6.129356080479549e-11
[-2.6875 -5.25  ]
-1 1.0115047418711676 5.270265888755516e-11
[-2.6875 -4.625 ]
-1 1.0008702989014906 5.306090409146694e-11
[-2.6875 -4.    ]
3937 0.0002932794012297095 3.411923390472684e-07
[-2.6875 -3.375 ]
3681 5.8403258277929774e-05 1.7187493838109818e-06
[-2.6875 -2.75  ]
2604 1.012427112227171e-05 9.912131377080882e-06
[-2.6875 -2.125 ]
2855 1.6052691211324852e-05 6.324432761515932e-06
[-2.6875 -1.5   ]
3478 0.00017115929424612184 5.901450571982897e-07
[-2.  -6.5]
-1 1.0040253690720715 2.6167606982977296e-12
[-2.    -5.875]
-1 1.0032273981557358 2.3352559702402662e-12
[-2.   -5.25]
-1 1.0060595489943398 1.0974989210872656e-12
[-2.    -4.625]
-1 1.001418119987236 1.0791395082965186e-12
[-2. -4.]
4277 0.0011589288607832251 8.655713815163773e-08
[-2.    -3.375]
3775 0.0003808965831836898 2.657351064981391e-07
[-2.   -2.75]
3148 7.834479864806533e-05 1.2773469597956814e-06
[-2.    -2.125]
-1 1.0008564664043758 1.789964616690148e-15
[-2.  -1.5]
3091 0.00030635218203736557 3.270999019032395e-07
Out[78]:
Text(0,0.5,'Normalized CO2 concentration')
In [79]:
"""
An alternative way of defining rates.
"""
new_rates = []
for experiment in experiments:
    #print experiment
    data = get_data(experiment)
    if data is None:
        print("No data for experiment {}".format(experiment))
        new_rates.append(np.nan)
        continue
        
    times = np.array(data.index)
    normalized = data[['CO2']].values[:,0] / highest_co2
    try:
        i = (np.nonzero((np.log10(times)+np.log10(normalized)) > -10))[0][0]
        time = times[i]
    except IndexError:
        time = 1.0
    new_rates.append(1./time)
new_log_rates = np.log(np.array(new_rates))
print new_log_rates
/Users/emilymazeau/anaconda2/envs/rmg_env/lib/python2.7/site-packages/ipykernel_launcher.py:16: RuntimeWarning: divide by zero encountered in log10
  app.launch_new_instance()
[ 7.23490282 10.65292844 13.33589734 13.60724387 13.83847104 15.03011758
 15.2180784  15.17183557 14.54754722  8.52871232 11.9108166  13.1707839
 13.05487007 13.6379215  14.77510029 14.83612262 14.77749621 13.51158467
  8.9844623  11.8679851  12.15393473 12.01927895 13.15028036 14.24872245
 14.26278295 14.05028323 12.25260654  8.62090753 10.84190583 10.78751737
 10.59966253 12.55895288 13.68966961 13.64867148 12.8888052  11.37110234
  8.53099856  9.22987783  8.99538106  8.84893998 11.61883249 12.80855854
 12.6303849  12.22419798 10.77871488  6.76287613  6.79575512  6.47605602
  6.25451123  9.91923565 11.32097702 11.76789357 11.85012325 10.01368994
  3.05162348  3.00267792  2.63830476  2.39800836  8.42528497  9.97402873
 11.62145434 11.44074855  9.37023141  0.          0.          0.
  0.          8.13438482  9.74813888 11.50057494 11.03963405  8.67291589
  0.          0.          0.          0.          6.7602591   7.87298265
  9.45439098  0.          8.0907752 ]
In [80]:
new_rate_grid = np.reshape(new_log_rates, (grid_size,grid_size))
plt.imshow(new_rate_grid.T, interpolation='none', origin='lower', extent=extent2, aspect='equal')
for metal, coords in abildpedersen_energies.iteritems():
    color = {'Ag':'w','Au':'w',}.get(metal,'k')
    plt.plot(coords[0], coords[1], 'o'+color)
    plt.text(coords[0], coords[1], metal, color=color)
plt.xlim(carbon_range)
plt.ylim(oxygen_range)
plt.xlabel('$\Delta E^C$ (eV)')
plt.ylabel('$\Delta E^O$ (eV)')
Out[80]:
Text(0,0.5,'$\\Delta E^O$ (eV)')
In [81]:
new_rate_grid = np.reshape(new_log_rates, (grid_size,grid_size))
plt.imshow(new_rate_grid.T, interpolation='spline16', origin='lower', extent=extent2, aspect='equal')
for metal, coords in abildpedersen_energies.iteritems():
    color = {'Ag':'w','Au':'w',}.get(metal,'k')
    plt.plot(coords[0], coords[1], 'o'+color)
    plt.text(coords[0], coords[1], metal, color=color)
plt.xlim(carbon_range)
plt.ylim(oxygen_range)
plt.xlabel('$\Delta E^C$ (eV)')
plt.ylabel('$\Delta E^O$ (eV)')
Out[81]:
Text(0,0.5,'$\\Delta E^O$ (eV)')

Count the reactions¶

In [82]:
def get_reaction_count(experiment):
    d = directory(*experiment)
    f = os.path.join(d,'chemkin','chem_annotated.inp')
    with open(f) as chemkin:
        r = re.compile('\! Reaction index: Chemkin \#(\d+)')
        for line in chemkin:
            m = r.match(line)
            if m:
                count = int(m.group(1))
    return count
    
get_reaction_count(experiment)
Out[82]:
373
In [83]:
i=35
print experiments[i]
print get_reaction_count(experiments[i])
[-5.4375 -1.5   ]
536
In [84]:
reaction_counts = map(get_reaction_count, experiments)
reaction_counts_grid = np.log10(np.reshape(reaction_counts, (grid_size,grid_size)))

plt.imshow(reaction_counts_grid.T, interpolation='none', origin='lower', extent=extent2, aspect='equal')
for metal, coords in abildpedersen_energies.iteritems():
    color = {'Ag':'w','Au':'w','Cu':'w'}.get(metal,'k')
    color='k'
    plt.plot(coords[0], coords[1], 'o'+color)
    plt.text(coords[0], coords[1], metal, color=color)
plt.xlim(carbon_range2)
plt.ylim(oxygen_range2)
plt.xlabel('$\Delta E^C$ (eV)')
plt.ylabel('$\Delta E^O$ (eV)')

for e,n in zip(experiments,reaction_counts):
    plt.text(e[0],e[1],n,color='w',ha='center', va='center')
#plt.colorbar()
In [85]:
# A linear one, just to check it looks the same
reaction_counts_grid = np.reshape(reaction_counts, (grid_size,grid_size))
ax = sns.heatmap(reaction_counts_grid.T[::-1,:], annot=True, fmt='d', square=True)
In [ ]: